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Abstract 

Thermal properties of low-density neutron matter are investigated by determinantal quantum 
Monte Carlo lattice calculations on 3+1 dimensional cubic lattices. Nuclear effective field the- 
ory (EFT) is applied using the pionless single- and two-parameter neutron- neutron interactions, 
determined from the ^Sq scattering length and effective range. The determination of the interac- 
tions and the calculations of neutron matter are carried out consistently by applying EFT power 
counting rules. The thermodynamic limit is taken by the method of finite-size scaling, and the 
continuum limit is examined in the vanishing lattice filling limit. The ^Sq pairing gap at T ~ 
is computed directly from the off-diagonal long-range order of the spin pair-pair correlation func- 
tion and is found to be approximately 30% smaller than BCS calculations with the conventional 
nucleon-nucleon potentials. The critical temperature Tc of the normal-to-superfiuid phase tran- 
sition and the pairing temperature scale T* are determined, and the temperature-density phase 
diagram is constructed. The physics of low-density neutron matter is clearly identified as being a 
BCS-Bose-Einstein condensation crossover. 
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I. INTRODUCTION 



Neutron matter is of great interest in nuclear physics as a quantum many-body system. 
The ^5*0 nucleon-nucleon (NN) interaction is strongly attractive, dominating the physics of 
neutron matter. The interaction yields the negative (in our convention) scattering length 
flo of an unnaturally large magnitude (~ 20 fm), with the effective range tq of a moderate 
(natural) size of about twice the pion wavelength (^ 2.8 fm). The value of implies that the 
strongly attractive interaction nearly forms a bound state. By this pairing, neutron matter 
is a strongly interacting many-body system, which must be treated nonperturbatively 

The strong neutron pairing generates a pairing gap that creates superfluidity in neutron 
matter. Superfluidity in neutron matter is of astronomical interest because of the close 
relation to the internal structure and thermal evolution of neutron stars 5*0 and 

^P2-^F2 superfluidity are believed to be realized in the inner crust and in the core region of 
neutron stars, respectively, and to contribute to the thermodynamic and dynamic properties 
of the stars. 

Neutron pairing is also considered important for understanding the structure of neutron- 
rich unstable nuclei. Neutron-neutron correlations are expected to be a crucial ingredient in 
the weakly bound, surface structure near the neutron drip line: and for the surface structure, 
neutron pairing in neutron matter must be well understood 

Investigations over many years have provided much understanding of the physics of ther- 
modynamic properties of neutron matter j3, ^,0], but reliable quantitative information of 
the thermal properties has not been fully available . For example, the ^5*0 pairing gap 
at zero temperature A had been firmly established in the BCS approximation, as evident in 
the fact that various conventional A^A^ potentials have provided nearly the same Ag^g as a 



function of neutron matter density 10|, lUl]- Many-body calculations beyond the BCS mean- 



field approximation, however, have yielded A of various magnitudes, generally smaller than 
the BCS values, some even by a factor of 2 or more. Quantum Monte Carlo calculations, 
based on a nonperturbative approach, have also been used on the A determination. The 
Green's function Monte Carlo (GFMC) method, quite successful in treating the ground-state 
properties of finite nuclei by the use of the conventional A^A^ potentials |12|, has yielded A 
in the low-density region {kp < 0.6 fm~^), smaller than Abcs fisl - 1^ but not as small as 
those obtained by some of the many-body calculations. Another method closely related to 
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GFMC, the auxiliary field diffusion Monte Carlo (AFDMC) method, which is also applied 
to finite nuclei [15], has given A quite close to Abcs 0, 12\ and significantly larger than 
the GFMC A. We present a more detailed comparison of these works, including ours, in 
Sec. Vll B. 

In this paper, we report a quantum Monte Carlo calculation of A and thermal properties 
of neutron matter using a method different from the GFMC and the AFDMC methods. The 
difference is that ours is based on the standard finite-temperature, grand canonical formu- 
lation, while the GFMC and AFDMC methods are based on essentially zero-temperature 
formulations, performed for the ground or specific excited states with a pre-fixed neutron 
number. Our calculation may be viewed, in a sense, as a nonrelativistic hadronic version 
of lattice QCD calculations, but it involves different aspects such as those associated with 
the large numbers of fermions on the lattice 1^. We use a Hamiltonian formulation differ- 
ent from the Lagrangian formulation commonly used in the lattice QCD calculations. Our 
brmulation is not new, as it has been applied in condensed matter physics for many years 



IS 



2m\ and has been also applied in nuclear physics j2]|. This work is an extension of the 



latter. 

We also use a new ingredient, the A^A^ interaction based on effective field theory (EFT) 



221 . |23[ |. in place of the conventional A^A^ potentials. It is desirable to include pions [24] 



in the EFT interaction as dynamical degrees of freedom, representing chiral symmetry and 
its breaking. Our objective is twofold: (1) to apply the A^A^ EFT interaction to the many- 
nucleon system of neutron matter by properly applying EFT counting rules, and (2) to 
determine reliably the thermal properties of neutron matter and their key quantities, such 
as A. In the first attempt for achieving this objective, we have chosen a pionless A^A^ EFT 
potential with two parameters. The major consequence of this choice is that application of 
our calculation is limited to the low-density region, kp < 0.6 fm^^. Even with this potential, 
our work has become a relatively large-scale computation, especially because we take the 
thermodynamic limit and examine the continuum limit. Note that field theoretical aspects 
of the general approach of this work were discussed a few years ago 25 ] . 

Because the pairing in neutron matter is strong, neutron matter should be treated as a 
strongly correlated fermionic system in the state of BCS-Bose-Einstein condensation (BEG) 
crossover, which has been receiving much attention in recent years [26]. Traditionally the 
pairing in neutron matter has been discussed in the framework of the BCS approximation 
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271], but the pairing is too strong for a BCS treatment. The pairing strength is characterized 



by l/lkpao) and corresporids to the BCS hmit with l/lkpCLo) —>■ — oo and to the EEC hmit 
with l/{kpaQ) +00 [28]. The range of l/lkpao) in the low-density region investigated 
in this work is well in the middle of the two limits, —0.8 < l/lkpao) < —0.1, and the 
magnitude of l/lkpao) becomes smaller for a higher density. We elaborate on the issue of 
crossover in Sec. VII A. 

The limit l/lkpcio) corresponds to the unitary limit, to which much attention has 
been paid lately in the fields of atomic and condensed-matter physics. A fermion pair in the 
unitary limit forms a zero-energy bound state, thereby yielding a scattering length infinitely 
long, associated with no classical scale and expected to have a universal feature. Our single- 
parameter EFT description of low-density neutron matter is close to the unitary limit (rather 
than to the BCS limit), and we will discuss the relation between the two in an accompanying 



paper [29|. We emphasize, however, that the close similarity of the two is restricted to the 
low-density region of neutron matter {kp < 0.3 fm~^), because additional EFT parameters 
and the pionic contributions needed for the description of the denser region introduce new 
length scales and make the physics more complicated than that of the unitary limit. 

The outline of this paper is as follows. After the Introduction of Sec. I, the basic setup of 
our calculation is described in Sec. II. In Sec. Ill, we present how we determine the physical 
quantities of interest in this work, and in Sec. IV, we show how we carry out their numerical 
calculation by taking the thermodynamic and continuum limits. In Sec. V, we discuss how 
the single- and two-parameter calculations are matched. The summary results are shown in 
Sec. VI, and discussions of the key points in this work are given in Sec. VII. A summary 
of our work is found in Sec. VIII. We include, in Appendix A, a relevant, short discussion 
on how the two A^A^ potential parameters are determined by satisfying EFT counting rules; 
in Appendix B, a comparison of the physical sizes of a neutron (Cooper) pair and the 
computational lattices; and, in Appendix C, somewhat detailed technical aspects of our 
Monte Carlo calculation. 
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II. BASIC SETUP 



A. NN EFT Hamiltonian 

The nuclear EFT Lagrangian is constructed by including all possible terms allowed by 
symmetries of the underlying theory of QCD 30|]. The A^A^ potential from the EFT La- 



grangian is written in the momentum expansion form 

Vip',p) = co(A) + C2(A)(p2 + p'') + ■■■- 2c2(A)p ■ p' + • • • , (1) 

where p and p' are the A^A^ center-of-mass momenta, and A is the regularization scale. The 
terms not explicitly shown in Eq. ([1]) include those in which pions are treated as a dynamical 



degree of freedom 3l|. For the momentum below the pion mass scale, we may neglect the 
explicit dynamics of chiral symmetry and its breaking by truncating Eq. ([T]) and including in 
Co and C2 the consequences of the dynamics. In this work, we use this pionless S'-wave A^A^ 
potential with the first two terms in Eq. ([1]). Generally an EFT potential is regarded as an 
expansion in terms of p/Q and p'/ Q with Q setting the momentum scale of the expansion. 
In our pionless potential, we have Q > {fn^^, the pion mass). Note that the potential 
consists of the central and spin-dependent parts, as Cc + cr ■ cr'co-, with a ■ a' = —3 for the ^Sq 
state (and = +1 for the ^Si state, not considered in this work). We also neglect in this work 
the P-wave interaction term starting with the p ■ p' and the relativistic effects appearing in 



0{p^/M^) |32|. 



Regularization is required for the application of Eq. ([T]). On a cubic lattice, the lattice 
spacing a serves as the regularization scale A, approximately as 

A ~ -. (2) 
a 

A should generally be set large, at least larger than the momentum p, 

A > (3) 



or better set 



A > Q, (4) 



corresponding to a < 4.5 fm for Q ~ 32, l33|. When the two-nucleon interaction is 



applied to a many-nucleon system of finite density, an additional constraint is imposed on 
the value of a, as discussed in Sec. II C. 



On the lattice, the Hamiltonian for our potential takes the discretized form [34 1 



1 

+ 



a3 



Co{a) ^C2(a) 

{l,j)(7<j' 



where t = l/(2Ma^), the hopping parameter (M is the neutron mass), and denotes a 
restriction on the sum to all neighboring pairs. and Qo- are the creation and annihilation 
operators of the neutron, with a =t, i, respectively, at the ith site. 

The neutron-neutron interaction parameters, Co{a) and 02(0), are determined from the 
neutron- neutron scattering phase shift, using the ^Sq effective range expansion (ERE), 

pcot 6o{p) = -— + Irop^ - Pry + 0{p^), (6) 

where P is the shape parameter. By dividing both sides by Q, we find Eq. ([6]) is an 
expansion in terms of the dimensionless quantity p'^/Q'^. For Q ^ ttIt^, the coefficients of 
the expansion roQ/2 and P(roQ)^ are of the natural size 0(1), while the first coefficient 
is unnaturally small, |l/aoQ| ^ 1. Phenomenologically the sum of the first two terms in 
Eq. ([H]) agrees well with the phase shift up to the center-of-mass momentum of nearly the 



pion mass ^0.7 fm ^, or about 40 MeV of the laboratory kinetic energy 35| (see also 



Ref. [36|]). This assures us that Co(a) and C2(a) are safely determined from ao and tq for a 



chosen value of a 



These interaction parameters are determined by consistently applying EFT power count- 
ing rules in a way different from a mere phenomenological fitting, as briefly discussed in 
Appendix A. Because this determination is one of the crucial steps in this work, let us note 
its key point here: C2(a) and the contributions of the same order must he treated perturbatively 
by neglecting the (9([c2(a)]^) contributions, so that 0{p'^/Q'^) contributions are consistently 
neglected. Furthermore, to be consistent, C2(a) and the contributions of the same order must 
also be treated perturbatively in the neutron matter calculations. In the next subsection, we 
discuss how this treatment is formulated for the neutron matter calculation. 

In this work, we carry out the neutron matter calculation using Eq. in two different 
ways: the leading-order (LO) calculation, in which the 02(0) contribution and the contri- 
butions of the same order are neglected, and the next-to-leading-order (NLO) calculation, 
in which they are included. The LO and NLO calculations are expected to yield somewhat 



different physics, because Eq. ([5]) is the Hamiltonian of the attractive Hubbard model for 
the LO calculation, and it is the Hamiltonian of an extended attractive Hubbard model for 
the NLO calculation 3J]. With the neglect of the LO calculation involves the 



neutrons of low momenta and should be applicable to a low-density region of neutron matter 
without the perturbative treatment. 

An important issue in this work is the density at which the LO and NLO results should be 
matched. The ERE of Eq. ([6]) suggests that the center-of-mass momentum of an interacting 
neutron pair is less than -^2/(1001^0) ~ 0.20 fm^^ at the matching density. As a rough 
estimate, it may be feasible ^o identify the Fermi momentum kp as this momentum and to 
estimate the density from it [1] , but for a rigorous matching, the LO and NLO neutron matter 
calculations should be carried out for some common densities and their results compared. 
As it is desirable to avoid excess computer time, we use in this work the following procedure: 
we carry out the LO and NLO calculations at the common density of kp = 0.3041 fm~^, 
where we expect the two results will certainly differ, and then perform similar calculations 
by lowering the density so as to identify the density that yields the same LO and NLO results 
(within the statistical uncertainties). The matching using this procedure is elaborated in 
Sec. V. 



B. Determinantal quantum Monte Carlo computation 



We follow a lattice Hamiltonian formulation, somewhat different from the Lagrangian 



formulation usually used in lattice QCD 



Instead of using the representation in terms 



of coherent-state Grassmann variables, we use the number representation, working with the 



lattice Fock space {n\ using the creation and anni 



21 



illation operators of the neutrons. Our 



39l | and is commonly used in condensed- 



treatment is the same as that used in Refs. 

matter physics [lol, [2^ under the determinantal quantum Monte Carlo (DQMC) method. 

We carry out neutron matter calculations using the Hamiltonian of Eq. ([5]) in the method 
of grand canonical ensemble. The Monte Carlo computation is carried out for various values 
of the chemical potential /i, and the /i dependence is converted to the density dependence 
by determining the densities by the average over i, a of {c\^Cio) for various values of /i. 

For many-nucleon systems, the Hamiltonian ^ should also include three-nucleon inter- 
actions. By EFT power counting rules, the interactions are to be treated generally as the 
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LO order in the pionless case, and they play a significant role when a three-nucleon bound 
state such as the triton can be formed 40|. In neutron matter, however, the three-neutron 
system has no bound state, and the three-neutron interactions appear at a higher order 
because the Fermi statistics prohibit the LO diagram of three neutron from being at the 
same spatial point with the momentum-independent vertex. As the interactions would also 
affect the two-neutron pairing indirectly, we expect that the interactions would affect the 
observables of our interest relatively weakly and defer the issue to a future investigation by 
neglecting them in this work. 

We write the partition function as 



Z{T,fi) = {n\Umn), 



(7) 



where U {(3) is the (imaginary time) evolution operator, and the trace implied in Eq. ([7]) is 
over all possible nucleon configurations on the lattice (n|. Using the Trotter- Suzuki approx- 
imation, we express U{I3) as 

Nt 



U{(3) = Texp 



Tt = l 



Tn^Lif/(A/3) 



(8) 



by the temporal discretization (3 = AjSNt, with A^^ being the number of time slices. In 
Eq. ([8]), -ff is the two-parameter NLO Hamiltonian of Eq. ([5]), and i is actually an integer 
vector specifying the location of a site with its component ranging as 

The Tt dependence of H and U{A/3) is solely through and c, as seen from Eq. ([5]). The 
last expression in Eq. (IHl) is thus a product of U{A(3) operators, each having the same form 
and depending on implicitly. 

To cast Z{T,fi) in a form amenable to Monte Carlo computation of the fermion integra- 
tion, we express the two-nucleon interaction of ^ in a single-nucleon interaction form by 
applying the Hubbard-Stratonovich transformation 



,+Ah'f 



— / dxiC 

71 



(9) 



for a constant A with Re(y4) > 0. Here, Xi is an auxiliary scalar field at the ith site, and fii 
is the density operator defined as fii = hi^ + hi\^ ( hi^ = c\^Ci^, the number operator with 
the spin a at the ith site). H is divided into two parts. 



H = 



(10) 



8 



where 



H' = — rCofa) fijn, H - 



{id) 



6 



Aco(a) -C2(a 



5Z ~ ^0 



(11) 



Here, Co(a) is expressed as a sum of the LO part c"Q\ci) and the NLO part Aco(a), which 
are defined in Eqs. flA3l) and flA2p . respectively, with A = n/a. 

We introduce Ho{x), the LO single-nucleon Hamihonian interacting with the external 
scalar field x = {Xi}? 

^o(x) = + J^X.n,. (12) 



In terms of Hq{x), U{A(3) is written as 
U{AP) = [ d[x\exp 



2a3 



exp <^ -A/5 



(13) 



where the measure is defined as d[x] = dxidx2 ■ ■ ■ with a constant factor generated by the 
Hubbard-Stratonovich transformation. We emphasize that H' is defined to be of the NLO 
and is treated perturbatively in the second step of Eq. f[T^ . 
We thus obtain 

= J d[x]G{x){n\U^mn), (14) 

where 



(15) 



Note that the time-ordering (sequential) integration over [x] is understood in the last ex- 
pression of Eq. f|T^ . In accordance with Eq. f|T3l) . the factor (1 — Aj3H') in U^{Aj3) of 
Eq. (fT5l) is to be evaluated by the use of the nucleon lattice configuration resulting from the 
exp{—Aj3[HQ{x) —fJ'J2i operation at r^, and that this procedure is repeated successively 
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from Tt = to Nf. This step is vital in the computation for the perturbative treatment of 
H'. We make a technically important note: because of the perturbative treatment of H', 
the number of the auxiliary fields for the NLO calculation remains as N^Nt, the same as 
for the LO calculation. If H' were not treated perturbatively, 4 x NgNt more of {xi} would 
have been needed owing to the derivative interactions, and the Monte Carlo computation 
would have required more time by nearly an order of magnitude. 

The trace of the single-particle evolution operator Uxif^) is expressed in terms of the 
single-particle matrix representation of the operator, U^i^P), as 



19 



20 



21 



39| 



(nlU^mn) = det [1 + C/^(/3)] ^ ^W- 



(16) 



The expectation value of the (static) operator (9(c^, c) at T = 1//? is then obtained from 

1 



(0(ct,c)) 



d[x]G{x){n\0{c\c)UM\^) 



where (C(x)) is 



z(r,/i) 

J d[x]Gixmx))ax) 

" fdixMxnx) ' 



(17) 



19 



20 



21 



{n\UM\n) 

and can be evaluated in terms of U^i^P) using Eq. (fT6|) . as shown in Refs. 

Equation f|T71) is now amenable to a Monte Carlo integration by treating \G{x) 
\G{x)^{x) \ as a weight. Our Monte Carlo computation is the same as that used in Ref. 



(18) 



39|. 



or 



supplemented by a matrix-decomposition stabilized method for low-temperature computa- 
tions 

Before closing this subsection, we make a relevant comment. In the procedure just de- 
scribed, we reduced the original Hamiltonian H of Eq. ([5]) to the single-nucleon Hamiltonian 
Hq (with H') of Eq. (fT2l) in terms of the density operators {ui}, as in Eq. (ITOl) . The choice 
of the density operators in this step may seem natural, but it is not required for the re- 
duction to an effective single-nucleon Hamiltonian because of the arbitrariness in the path 
integral formulation. In fact, we can choose a combination of pairing operators and density 
operators, leading to a Hartree-Fock-Bogoliubov (HFB) type Hamiltonian jsQ, 41|. 
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C. Lattice spacing toward thermodynamic and continuum limits 



Neutron matter is a strongly correlated fermion system. On a three-dimensional cubic 
lattice, the correlation length resulting from the simulation, ^, satisfies 

a « e < (19) 

1/3 

where L = aNg is the physical dimension of the cubic lattice. ^ is the length scale in which 
the collective state is realized in the simulation and is different from the size of a neutron 
pair (a Cooper pair) in the state ^cp- Note that, confusingly, has often been referred 
to terms similar to ^. In Appendix B, we compare the physical sizes of the neutron pair 
simulated and the lattice spaces used. 

To obtain a physically meaningful result, we seek for and for the expectation values of 
other quantities, in the continuum limit a —>■ and in the thermodynamic limit L oo. 
The clear procedure for achieving both limits is to do the former with L fixed (for obtaining 
results insensitive to the lattice structure), and then to do the latter (using finite volume 
corrections), as is usually done in lattice QCD calculations 



In our calculation of the many neutron system, each meaningful configuration must consist 
of neutrons fewer than A^^, so that the calculation properly describes the interacting neutron 
system in free space. This requirement is crucial in general for the simulation of a system of 
many fermions, and we find that the requirement complicates the straightforward approach 
of achieving the above two limits. Note that lattice QCD calculations have not yet dealt 
with cases of such high baryon-density states. 

Let us elaborate on this requirement. Consider setting up a classical lattice configuration. 
When Nf neutrons are placed on a lattice of volume a^Ng, the neutron density p is 



Nf n 



which defines the lattice- filling fraction (or more descriptively, the site-occupation fraction), 
n = Nf/Ng. n denotes the fraction of the lattice sites occupied by the neutrons. Note that 
the complete filling of the lattice occurs with n = 2 owing to the spin degree of freedom. 
Classically, n can simply be chosen, while in our quantum-mechanical, grand canonical 
calculation, it is determined from Xli cr(^L^«o-); "which is computed for a fixed value of a and 
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Mathematically, for a finite nucleon density, Eq. (120!) implies 

(21) 

as a ^ 0. Physically, these limits simulate the free-space environment, because the smaller 
n is, the more vacant sites are available, allowing more feasible excitations to be realized. To 
determine thermal quantities as a function of neutron matter density, we consider achieving 
the limits to be vital and take the limit of Eq. ( 1211) as the continuum limit. Note that this 
procedure is similar to, but different from, the one recently proposed for the unitary limit 



problem 42|, in that we keep the density p finite as we approach the continuum limit, but 



the /c/r — *^ limit is taken in Ref. 



42| 



Once we decide to take Eq. (12T]) as the continuum limit, we have to use different values 
of a for different densities to satisfy the regularization scale requirement, Eqs. ^ and (jl]), 
of the EFT. The procedure becomes complicated in order to satisfy all these requirements, 
but at the same time it has to be durable in practice. We have decided to use the following 
procedure. First, we choose an appropriate n value that is small enough yet reasonably 
durable. Second, for this n, we choose a set of the representative nucleon densities for 
the computation and a set of appropriate a values for them. We call the set the standard 
parameter set, and we list them in Sec. II D. Third, after we complete the computation for 
the standard set, we perform the computation by varying the lattice size, so as to take the 
thermodynamic limit. Fourth, we vary n to examine the continuum limit as n ^ 0. 

In the rest of this subsection, we discuss the first step, how we choose n for the standard 
set. As an estimate, take the Fermi gas model. In terms of the Fermi momentum kp, n for 
neutron matter is written as 

n = {akpf/i^n'^) 0. (22) 

To keep n independent of a for various densities, we should have a oc l/kp- Note that the 
excitation energies of the neutron matter of interest are about an order of magnitude less 
than the Fermi energy, as seen in Sec. V, and are safely ignored in this estimate. 

The smallness of n is achieved by making a small, or A large. If we take Eq. ([3]), Eq. ([2]) 
with p kp yields 

TT > akp. (23) 
Equations ( l22i) and (123|) yield a rather loose estimate of n < 1. We can obtain a more realistic 
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limit from the observation that the lattice discretization amounts to the replacement 

2 ^ 

- - E [1 - = P' + ^(«'P')' (24) 

i=l 

for example, in the neutron propagator. This observation suggests that the left-hand side 
of Eq. (1231) is more like unity instead of vr, and we obtain the inequality 

(37r2)-i/3 > n. (25) 

This choice of n does not require a large A to satisfy Eq. (jlj), but it does for Eq. 

The preceding consideration leads us to set n = 1/4 (or 1/8 of the full filling of the 
lattice), as a practical compromise. Other parameters also need to be chosen. In the 
following subsection, we discuss how they are chosen and list all parameter values in the 
standard parameter set. 

D. Standard parameter set 

The standard set of the potential parameters is shown in Table HI We choose the set by 
the following steps. First, we choose the values of the Fermi momentum kp, representing 
the neutron matter density, as shown in the first and second columns. Second, the values 
of a are determined from {akF)^/{3'n-'^) = n = 1/4 and at the same time by ensuring that 
the a values provide reasonable EFT regularization scales. Third, the values of Cq and C2 
are determined from oq and tq using the a values in Eqs. (]A2p and (1A3I) with A = vr/a. The 
Monte Carlo calculations are carried out using the cq and C2 values by tuning the chemical 
potential fi, so that the resultant neutron matter densities by the Monte Carlo computation 
are the p values listed in the third column in the unit of the normal nuclear density po = 0.16 
fm~^. We emphasize that these p values are expressed in terms of the kp values in the first 
and second column exactly as 

p = 4/(37r2). (26) 

Throughout this work, we use kp defined through Eq. ( 126|) for specifying the quantum- 
mechanically computed density, p, of neutron matter as the interacting fermion system. 
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TABLE I: Standard parameter set. 



kp 


(MeV) k 


F (fm-1) 


P {po) 




a (fm) co/(a^f) 


C2/{aH) 


LO 


15 


0.07602 


9 X 10" 


5 


25.64 -5.308 


- 


LO 


30 


0.1520 


7 X 10- 


4 


12.82 -6.354 




LO 


60 


0.3041 


6 X 10- 


3 


6.409 -7.049 




NLO 


60 


0.3041 


6 X 10- 


3 


6.409 -9.646 


0.3684 


NLO 


90 


0.4561 


2 X 10- 


2 


4.273 -11.074 


0.5139 


NLO 


120 


0.6081 


5 X 10- 


2 


3.205 -12.343 


0.6555 



III. DETERMINATION OF A, T^, AND T* FROM THE PAIRING CORRELA- 
TION FUNCTION 

In this work, we focus on the determination of three quantities: the pairing gap at 
T ^ 0, A; the critical temperature Tc of the normal-to-superfluid phase transition; and the 
pairing temperature scale T*. The latter two will be used to obtain the density-temperature 
phase diagram, and all quantities will be calculated from correlation functions, the first two 
from the pair-pair correlation function and the third from the magnetic susceptibility (the 
spin-spin correlation) . 



A. Pairing gap A 

A is determined directly from the off-diagonal long-range order (ODLRO) of the spin 
pair-pair correlation function Pg 43|, 



i 

= ^ E ^^^^ - ^^■^)' ' (27) 

i, j=i+R 

where Aj = QiQi is the two-neutron spin-pairing operator at the iih. site, and R is the 
separation of the neutron pairs in the lattice spacing unit and has no dimension. Note that 
Gij = Glj = (ci^c^j^) for a =|, | in the attractive Hubbard model. Ps{R) decays rapidly 
in i? ^ 1 or 2 and takes a diminishing asymptotic value. When a long-range order exists 
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between neutron pairs, the asymptotic value is finite, that is, the signature of the ODLRO. 
Figure [1] illustrates this behavior. 

In Fig. [H Ps{R) is calculated for 14 values of T/t between 2.0 and 0.0625; but for clarity, 
only the selected values of T/t are shown. Note that the integer points of R = 1-4 arise 
from the lattice points in the side of the cubic, while the largest R = 4^/3 ~ 7 comes 
from the midpoint of the diagonal line in the cubic, which has the displacement vector 
(4,4,4). The values for i? > 3 are found to be quite close to each other at the lowest three 
temperatures, T/t = 0.25, 0.125, and 0.0625. The values at i? = 4 and 7 are averaged, 
yielding Ps(T ^ 0, R ^ 1). A is then determined from 

A = J^VP,(T^O,i?> 1). (28) 

Similar procedures are applied for different Ns and kp. The errors from the fit hereafter are 
estimated by a constrained least-squares method. 

"1 1 1 1 1 1 1 1 

T/t = 2 — I ■ 
T/t = 0.444 H-H — I ■ 
T/t = 0.25 ■ 
T/t = 0.125 I A - 

m 

— f-— g i-— : 

e s : 

Dm m ■ 
J i ^ ^ \ \ ^ 1 

012345678 

R 

FIG. 1: Spin pair-pair correlation function Pg as a function of the lattice separation R for the 
lattice size Ng = 8^ for the site-occupation fraction n = 0.25 at kp = 30 MeV. The DQMC results 
are shown with statistical uncertainties for T/t = 2.0, 0.444, 0.25, and 0.125 in the unit of hopping 
amplitude t = 0.126 MeV. The dashed line is the asymptotic value of Ps = 0.0244(6) extracted 
from the values for i? = 4 and 4^3 at T/t = 0.25, 0.125, and 0.0625 (not shown). 

Note that, as seen in Fig. [2] of the next subsection, the critical temperature is Tc/t = 
0.335(1), and the behavior of Ps{R) at Tc is similar to T/t = 0.25 in Fig. [H We caution the 

15 



0.06 
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0.02 



reader that the A thus determined is not our final value but is the value for Ns = 8^ and 
= 1/4 at fci? = 30 MeV. Using A's for various values of Ns and n at each kp, we determine 
A at the thermodynamic and continuum limits by the further analysis described in Sec. IV. 
The same caution is applied to the determination of and T* in the following subsections. 



B. Critical temperature Tc 



Tc of the normal-to-superfluid phase transition is determined from the spin pair-pair 



correlation sum 



44 



45 



46 



47| 



(29) 



Tc is extracted from the infiexion point of Ca(T). Figure [2] illustrates a typical case of Ca 
as a function of T/t, for kp = 30 MeV and Ns = 8^. In the figure, the infiexion point is at 
Tjt = 0.335(1), or = 0.0423(1) with t = 0.1261 MeV. The infiexion point is determined 
by an interpolation that fits the Monte Carlo data with an assumed function. 



Ca(T) = -Ci tanh[C2(T - T,)/t] + C3 



(30) 



where Ci = Ca(T = 0)/2, C2, and C3 are free constant parameters. 



C. Pairing temperature Scale T* 

As the temperature increases, the long-range order of the superfiuidity disappears at 
Tc. Above Tc, the spin pairing still remains, however, without generating the long-range 
order, and as the temperature increases further, the pairing eventually disappears. Though 
the process of the pairing disappearance is expected to be a continuous process, we may 
identify the temperature below which the pairing can be viewed as still strong. Following 



a practice in condensed-matter physics 



43], we denote the temperature as the pairing 



temperature scale T* and determine it from the temperature dependence of the Pauli spin 
susceptibility xp- When the (S'-wave singlet) spin pairing is weakened, the spectral weight 
of low-energy spin excitations is reduced, and the spin response weakens, xp is a good 
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FIG. 2: Spin pair-pair correlation sum Ca as a function of temperature T in the unit of the 
hopping parameter t. The Monte Carlo data with the statistical uncertainties are shown for the 
case of Ns = 8^ and kp = 30 MeV. The vertical dashed hue signifies Tc{Ns = 8^) = 0.335(l)t, 
corresponding to the inflexion point of the interpolated curve of C^{T/t), Eq. ([30]) with Ci = 
13.6 lb 1.0, C2 = 15.1 lb 2.9, and C3 = 16.2 ib 0.6. The interpolated curve is shown as the solid 
curve. 




quantity for studying this transition, since the xp of ^ fr^e fermion gas diverges as T — 0, 
while it vanishes for an interacting fermion g 3iS, clS illustrated in Fig. [3l 
Xp is given by 

id 

1 1 X 

^ 2^ ]\r ^ ^^'^ ~ ^^'^ ' ^^^^ 

where Si = c|^cr/^i/Cjv and cr is the Pauli vectorial matrices. T* is determined by 

identifying the maximum point of Xp as a function of T as discussed in the following. 

Figure m shows a typical case. For kp = 30 MeV and N = 8'^, we obtain T* /t = 0.5253(3) 
[T* = 0.06624(3) MeV with t = 0.1261 MeV]. The maximum point of the Monte Carlo data 
is determined through interpolation by use of a fitting function with a parameter Ci, 

Xp{T) = C^Texp{-T/T*). (32) 
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FIG. 3: Pauli spin susceptibility xp ^ a function of temperature T in the unit of hopping parameter 
t for Ns = 4'^. The sohd curve is the free fermion gas hmit (|co|/(a^t) 0) oixp{T), ~ n{l—n/2)/T 



(with the fining fraction n) 



48l |. In comparison to this, the cases for |co|/(a^t) = 2, 4, 6, 8, 10, and 



12 are shown in the increasing order of the interaction strength, from top to bottom. 
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FIG. 4: Pauh spin susceptibihty xp as a function of temperature T in the unit of hopping parameter 
t for Ns = 8^ and kp = 30 MeV. The vertical dashed fine signifies T*{Ns = 8^) = 0.5253(3)t, 
which is determined as the maximum point of xp{'^)i using the fitting function Eq. (j32|) with 
Ci = 0.258 ± 0.008 (shown as the sohd curve). 
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Note that though the definition of T* is somewhat subjective, T* thus defined approaches 
Tc at the BCS hmit, and T* signifies the pair-forming temperature at the BEC hmit as 
T* oc |co|/ [aH\n{\cQ\/{aHeF)f'^] Q, Q]- Here, the BCS and BEC hmits correspond to 
the weak and strong interaction hmits, or the small and large Co/(a'^t) limits, respectively. 



IV. A, Tc, AND T* AT THE THERMODYNAMIC AND CONTINUUM LIMITS 



A. Pairing gap A 



> 

CD 



1 

0.8 
0.6 
0.4 
0.2 h 



kp = 60 MeV 
kp = 30 MeV 
kF = 15MeV 



0.05 



0.1 



N, 
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FIG. 5: Lattice-size {Ns) dependence of A at LO for kp = 15, 30, and 60 MeV with n = 1/4. The 
Monte Carlo data shown with the statistical uncertainties are obtained for Ns = 4^, 6^, 8^, and 



10^. The dashed lines are the best fits by the use of linear functions of A'^. 



-1/2 



As the first step, we determine A at the thermodynamic limit. To carry out a definite 
analysis, we apply the BCS finite-size scaling exponent, A = 3/2 in A ~ = Ns as 
being independent of the density {49]. The exponent is obtained through A(T = 0, V^) oc 
Tc(V,) ~ = V7^^^ by combining the BCS resuh, A(r = 0) ^ 1.76Tc, and the direct 

relations between the finite-size scaling and critical exponents Q, [s^ . Note that while the 
usual best fit to all of our Monte Carlo data results in an essentially indefinite A, the 
jackknife method (often used in the lattice QCD data analysis [38]) yields A = 1.6 ± 0.3 
by assuming a linear dependence independent of the density. Apparently the value of 
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FIG. 6: Same as Fig. El but at NLO for kp = 60, 90, and 120 MeV. 



the exponent changes httle between the BCS weak-couphng region and the neutron-matter 
BCS-BEC crossover region. Figures [5] and [6] show the choice of A = 3/2 reasonable. 

In Figs. [5] and [6], the finite-size scaling of A is shown as a function of A^^ using A^^ = 4^, 
6'^, 8'^, and 10'^ data with n = 1/4. Fig. [5] is the finite-size scaling of A evaluated at LO for 
kp = 15, 30, and 60 MeV, while Fig. [6] is at NLO for kp = 60, 90, and 120 MeV. The data 
at LO are found to be best fit with a linear dependence on N^^^'^ = as 

A{Ns,kF = 15 MeV) = 0.0394(34) A^'^^/^ + 0.019152(20), 

A{Ns, kp = 30 MeV) = 0.096(35) N^^/^ + 0.1207(16), (33) 

A{Ns, kp = 60 MeV) = 0.74(24) N'^^^ + 0.581(13), , 

and those for NLO are 

A(iV„ kp = 60 MeV) = 0.423(67) N;^/^ + 0.4602(54), 
A(iV„ kp = 90 MeV) = 1.16(17) N;^^^ + 1.028(14), (34) 
A{Ns, kp = 120 MeV) = 3.91(75) N^^^^ + 1.565(42), 

where the last constant for each value of kp is A at the thermodynamic limit {Ng oo). 
The best-fit constants in Eqs. (1331) and (IMl) are determined using the jackknife method. 

As the second step, we determine A in the continuum limit using the above thermody- 
namic limit values. As discussed in Sec. II C, these values are obtained by using the standard 
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FIG. 7: Pairing gap A in the unit of the Fermi energy ep as a function of the filhng fraction n for 
kp = 60 MeV. The sohd circles show Monte Carlo data for Ng = 6^ at LO, with statistical uncer- 
tainties. The dashed line is the best fit by the use of a linear n^^^ dependence. The interception of 
the dashed line with the vertical axis corresponds to A at the continuum limit (n 0) for Ng = 6^ 
at kp = 60 MeV. 



parameter set, or with n = 1/4 (half of the quarter-filling), and are needed to extrapolate 
to n = to reach the continuum limit, a 0. For the extrapolation, we need to know how 
much A changes between n = 1/4 and n — > 0, or the ratio of A at the two values of n, Ra- 
In this work, we determine Ra solely using LO Monte Carlo data of the A^^ = 6'^ lattice for 
kp = 60 MeV. Dependence of Ra on A^^ aiid kp is weak both for LO and NLO, as discussed 
in Sec. IV D. 

Figure [7| shows the dependence of A on n for A^^ = 6^ at fci;' = 60 MeV. The data in 
the figure, shown with statistical uncertainties by solid circles, are for n = 1/16, 1/8, 3/16, 
1/4, 3/8, and 1/2. The EFT potential parameter co{a) is varied by the use of Eq. ( 1A3I) to 
accommodate the variation of a generated by the change of n. 

The n dependence of A is found to be relatively weak, and the jackknife analysis of the 
data yields 

A(n, Ns = 6^)/ep = -0.07(7) n^-^^^-^^ + 0.337(20). (35) 

While more data are desirable to reduce the uncertainty of the continuum limit, the constant 
term in Eq. (l35l) . some indirect information of the n exponent is available from the weak- 
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coupling BCS theory by the use of A oc Tc, and also from the analysis by Burovski et al. [42 1 
in a similar limit (but with kp —>■ a.s noted in Sec. 11 C) for their unitary limit calculation. 
Both suggest the n^^^ dependence, with which we find the best fit 

A(n, A^, = 6^)/eF = -0.044(16) n^^^ + 0.351(10). (36) 

For definiteness and because of lack of time, we use in our present analysis Eq. (136|) and show 
it as the dashed line in Fig. [71 Equations ( 135|) and ( l36l) yield the statistically consistent A at 
the continuum limit and suggest the systematic uncertainty by the use of the n^^^ dependence 
to be several percent. 

Equation fl36|) gives the ratio i?A 

_ A(n^0,7V, = 6^) _ 0-674(19) MeV _ 

- A{n = 0.25, N, = 6^) " 0.628(11) MeV " ^^^^^'^ ^'^^^ 

That is, the continuum- limit correction amounts to a 7% increase in the value of A. Ex- 
ploiting the weak dependence of on Ng and kp (elaborated in Sec. IV D), we apply the 
same -Ra to A at the thermodynamic limit in Eqs. ( !33l) and f l34l) . so as to obtain the final 
values of A at the thermodynamic and continuum limits. 



B. Critical temperature Tc and pairing temperature scale T* 



To obtain and T* at the thermodynamic and continuum limits, we carry out the same 
two steps as those done on A in the preceding subsection. Because Tc is at criticality, we will 
apply the universality argument for taking the thermodynamic limit. Monte Carlo data at 
n = 1/4 used for the finite-size scaling of Tc and T* are shown in Figs. [8] and [91 for A^^ = 4^, 
6'^, 8'^, and 10^ with statistical uncertainties. 

The exponent of the finite-size scaling and the critical exponents are known to be directly 
related at criticality 50|. Furthermore, because the three-dimensional (3D) XY model and 



our 3D Hubbard model are expected to belong to the same universality class 
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5l|, 



the exponents of finite-size scaling at criticality of both models are also expected to be the 



same [49|, [50|]. Accordingly, we have Tc{kF-,Ns) — Tc{kp,Ns — > oo 



A^. 



-1/(3^) 



with z/ = 2/3 of the XY model 5l|, IS^- Here, v denotes the exponent of, for example, the 



correlation length, as ~ (T — T^) . Note that in comparison, a mean-field approximation 



such as Ginzburg-Landau theory gives v = 1/2 [53|. With the linear A^^ dependence, we 
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find the best fits to the data for Tc at LO to be 

Tcikp = 15 MeV,N,) = 0.039(14) N'^^^ + 0.00700(94), 

T^ikp = 30 MeV, A^,) = 0.1839(11) iV^^/^ + 0.03420(11), (38) 

Tcikp = 60 MeV, AT,) = 0.6069(64) N;^/^ + 0.15770(30), 

which are shown in Fig. [HI and at NLO, 

Tcikp = 60 MeV,Ns) = 0.88(19) N;^/^ + 0.146(13), 
T^(kF = 90 MeV, Ns) = 1.14(33) N^^^"^ + 0.388(22), (39) 
T^ikp = 120 MeV,Ns) = 1.67(63) N^^^^ + 0.687(37), 

which are shown in Fig. [9l The last constant in each best fit in Eqs. (138!) and (139!) is at 
the thermodynamic hmit, Tc{kF,Ns oo). 

While T* is not at criticality, we find the finite-size scaling for T* to be similar to that of 
Tc. For example, the data of T*{Ns) yield the best-fit scaling power ~ jV-o-507±o.oo7 ^j^]^ ^^j^g 
jackknife method (for T^, ~ jy-o.53±o.o3^_ thus apply the same linear Ns ^^"^ dependence 
to T* as that for Tc. The best fits for T* at LO are found to be 

T*{kF = 15 MeV, AT,) = 0.1400(24) N;^/^ + 0.00707(24), 

T*{kF = 30 MeV, Ns) = 0.448(54) A^^^/^ + 0.0463(23), (40) 

T*{kF = 60 MeV, Ns) = 1.45(15) N^^^^ + 0.2618(99), 

and at NLO, 

T*{kF = 60 MeV,Ns) = 1.890(59) N^^^^ + 0.2575(35), 
T*{kF = 90 MeV, A^,) = 3.583(60) N^^^^ + 0.5876(61), (41) 
T*{kF = 120 MeV,Ns) = 3.70(12) N^^^^ + 1.4320(69), 

where the last constant in each equation gives T* at the thermodynamic limit. 

As to the continuum limit, in Fig. [10] we show the n dependence of Tc and T* at LO for 
A^s = 6'^ at fciT' = 60 MeV. The data with statistical uncertainties are shown by solid circles 
for n = 1/16, 1/8, 3/16, 1/4, 3/8, and 1/2. The exponent fit of (T*) shows ~ nO-3i±o i2 
(T* ~ ^o.43±o.io^_ observed for the similar limit of Tc [42|], they appear to be best fit by 
a linear n^^^ dependence, 

Tc{n, Ns = 6^)/eF = -0.165(23) n^^^ + 0.209(12), (42) 
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FIG. 8: Finite-size scaling of the critical temperature and the pairing temperature scale T*. 
The Monte Carlo data for Tc and T* with n = 1/4 at LO are shown for kp = 15, 30, and 60 MeV 
from bottom to top. The dotted lines are the best fits of Eqs. (f38l) and (HO]l . 
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FIG. 9: Same as Fig. [8l for kp = 60, 90, and 120 MeV, from bottom to top. The dotted lines are 
the best fits of Eqs. ([39]) and (liTI) . 



and 

T*(n, N, = 6^)/eF = -0.286(20) n^^^ + 0.367(12). (43) 
Note that the continuum limits of and T* in Eqs. ( H2|) and (H3|) are consistent with 
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FIG. 10: n dependence of the critical temperature Tc and pairing temperature scale T* at LO in 
the unit of the Fermi energy ep for Ns = 6^ at kp = QO MeV. The lines are the best fits to the Tc 
and T* data, Eqs. (^2!) and (H3]l . respectively. 

those determined by the exponent fits using the jackknife method within the statistical 
uncertainties [Tc = 0.223(41) and T* = 0.328(34)]. Contrary to the case of A, the n 
dependence of Tc and T* is rather strong. Equations (H2|) and (H3l) provide the needed ratios 
Rt^ and Rt*, which are used to obtain Tc and T* at the thermodynamic and continuum 
hmits, as in the case of A. 

C. Dependence of the continuum limit on Ng and kp 

The extrapolation to ri — > depends generally on A^^ and hp, but the dependence is 
expected to be weak because of the separation of local (ultraviolet) and global (infrared) 
properties for a sufficiently large A^^- 

For the A^^ dependence, we calculate, using the lattice sizes of A^^ = 4'^ and 8^, the ratios 
between n ^ and n = 0.25: -Ra, Rtc^ and -Rt*, both at LO and NLO. As summarized 
in Table HTl each ratio at hp = 60 MeV is consistent within the statistical uncertainties for 
Ns = 4^, 6^, and 8^ both at LO and NLO. Note that the second row for A^^ = 6'^ is obtained 
using data at n = 1/16, 1/8, 3/16, 1/4, 3/16, and 1/2, while the other rows use data at 
n = 1/16, 1/4, and 1/2. 
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Table IIIII also confirms the weak dependence on kp. Note that the third row uses data 
for n = 1/16, 1/8, 3/16, 1/4, 3/8, and 1/2, while the other rows use data at n = 1/16, 1/4, 
and 1/2. 



TABLE II: Dependence of the continuum limit on A'^; 





kp (MeV) N, 


Ra 


Rn 


Rt" 


LO 


60 


43 


L14(17) 


2.10(15) 


1.9(2) 


LO 


60 


63 


L07(5) 


1.96(13) 


1.94(9) 


LO 


60 




L12(8) 


1.86(10) 


2.0(1) 


NLO 


60 


43 


L08(6) 


2.08(8) 


2.1(2) 


NLO 


60 


63 


1.04(5) 


2.05(7) 


2.1(1) 


NLO 


60 


83 


1.05(4) 


2.08(37) 


2.0(1) 



TABLE III: Dependence of the continuum limit on kp. 





kp (MeV) Ns 


Ra 


Rt. 


Rt- 


LO 


15 


63 


1.09(3) 


1.96(9) 


2.4(9) 


LO 


30 


63 


1.08(8) 


1.98(6) 


2.4(4) 


LO 


60 


63 


1.07(5) 


1.96(13) 


1.94(9) 


NLO 


60 


63 


1.04(5) 


2.05(7) 


2.1(1) 


NLO 


90 


63 


1.11(10) 


2.12(20) 


2.0(1) 


NLO 


120 


63 


1.04(2) 


2.00(36) 


2.0(3) 



V. MATCHING LO AND NLO RESULTS 

Figures [TT] and [T2] display the LO and NLO A's as a function of kp and illustrate their 
matching in the region of kp = 0.15-0.30 fm~^: the A shown in Fig. [TT] is the result of 
the elaborate calculation described in Sees. Ill and IV, while the A shown in Fig. [12] is the 
result of a simpler calculation for 43 lattices with n = 1/4, including A at the density of 
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kp = 0.22805. The density dependences of the A's are quite close to each other in the two 
figures, demonstrating that a smooth transition from the LO A to the NLO A occurs in the 
density region of kp = 0.15-0.30 fm"\ Accordingly, we take the LO A for kp = 0.1520 fm"^ 
and the NLO A for kp = 0.3041 fm^^, as the final values. 
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FIG. 11: ^5o pairing gap, A, in the thermodynamic and continuum limits, resulting from the LO 
(solid circles) and NLO (open circles) calculations. The neutron density is denoted in terms of the 



Fermi momentum kp. The BCS calculation of Ref. 
including polarization effects of Ref 



54| (solid curve) and a higher order calculation 



55l | (dashed curve) are also shown for comparison. For a more 



detailed comparison, see Fig. [T7] in Sec. IVIIBi 



Figure [13] shows that also for and T*, smooth transitions take place between the LO 
and NLO values in the same density region as for A. We thus also take Tc and T* at 
kp = 0.1520 fm"^ as the LO and and T* at kp = 0.3041 fm"^ as the NLO. Note that 
the difference between the LO and NLO values of and T* in Fig. [12] is much smaller than 
that in the case of A. 



VI. RESULTS 



A. Pairing gap A 

Table IIVI lists our final values of A in the thermodynamic and continuum limits for 
low-density neutron matter. Table [IV] includes the ratio of A and the corresponding BCS 
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FIG. 12: pairing gap, A, for Ng = 4^ and n = 1/4, resulting from the LO (solid circles) and 
NLO (open circles) calculations. The neutron density is denoted in terms of the Fermi momentum 
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FIG. 13: Critical temperature Tc (circles) and the pairing temperature scale T* (squares) by the 
LO (solid symbols) and NLO (open symbols) calculations for Ng = 4^ and n = 1/4, shown as a 
function of the neutron matter density (represented by the Fermi momentum kp) The error bars 
are statistical uncertainties only. 
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TABLE IV: Our final values of the 5*0 pairing gap A in the thermodynamic and continuum limits, 
and the ratio of A and the BCS value Abcs- Uncertainties are statistical only. 



kp (MeV) 


P {po] 




A (MeV) 


A/Abcs 


15 


9 X 10" 


-5 


0.021(1) 


0.69(3) 


30 


7 X 10- 


-4 


0.13(1) 


0.67(4) 


60 


6 X 10" 


-3 


0.49(3) 


0.56(5) 


90 


2 X 10" 


-2 


1.10(7) 


0.68(4) 


120 


5 X 10" 


-2 


1.7(1) 


0.74(4) 



pairing gap, Abcs- Here, the Abcs's are taken from those tabulated in Ref. |5J] as the 
representative BCS values. As noted in Sec. VII B, there are only quite small differences 
among the Abcs calculated by the CD-Bonn, Nijmegen I, Nijmegen II, and Argonne V18 



A^A^ potentials [lo|,lll|- 



It is difficult to assess the systematic uncertainties involved in our calculation. In view 
of the probable uncertainties involved in taking the thermodynamic limit and especially 
the continuum limit, however, it would be fair to state that our calculation yields A to be 
approximately 30% less than the BCS values, perhaps with an additional systematic uncer- 
tainty of about ±10%. We thus consider finer variations of A inconclusive. For example, a 
close examination of Table HVl shows that the A/ Abcs ratio dips at around kp = 60 MeV. 
But this would require further study. 



B. Phase diagram of low-Density neutron matter 

Table |V] lists our final values of Tc and T* in the thermodynamic and continuum limits. 
It also shows their ratios and the ratios with the A of Table IIVI In Table IVj we observe 
that T* approaches Tc as the density decreases. That is, the pseudogap state (see below) 
diminishes as the density decreases. Furthermore, as the density decreases, the A/T^ ratio 
approaches the BCS value of about 1.76 [56], while A and themselves remain different 
from the BCS values. 

Tc and T* in Table IVl provide the temperature-density phase diagram as shown in Fig.[T4l 
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The figure illustrates the thermodynamic properties of low-density neutron matter. For 
example, at a Bxed den.ty k,, a. the tempetatute goe. dow^ttom the notmal phase, the 
pairing is gradually enhanced, forming the pseudogap phase [49] around and below T* . As 
the temperature goes down farther, the pairing gets stronger and eventually forms a long- 
range ordering at T^, thereby generating the second-order phase transition to the superfiuid 
phase. Note that the transition between the pseudogap phase and the normal phase is 
smooth. We must also note that the definition of T* is somewhat subjective. 

TABLE V: Our final values of and T*, and the relative magnitudes among them and A in Table 



kp (MeV) Tc (MeV) 


T* (MeV) A/Tc 


A/T* 


Tc/T* 


15 


0.014(3) 


0.014(1) 1.5(4) 


1.5(2) 


0.99(28) 


30 


0.067(5) 


0.091(9) 1.6(2) 


1.4(2) 


0.74(12) 


60 


0.29(5) 


0.45(5) 1.7(4) 


0.99(11) 


0.57(12) 


90 


0.76(9) 


1.1(1) 1.5(3) 


0.97(11) 


0.67(12) 


120 


1.4(2) 


2.8(1) 1.2(2) 


0.60(7) 


0.49(8) 



VII. DISCUSSIONS 

A. Nature of low-density neutron matter: BCS-BEC crossover 

To understand the nature of low-density neutron matter, we examine the dependence 
of Tc on the parameter cq by applying the LO calculation, since the physics throughout 
our low-density region is largely dictated by cq. Figure [15] illustrates the dependence in 
comparison to Tc in the weak-coupling (BCS) and strong-coupling (BEC) limits. 



T,(BCS) = — 7(36^2 -;,2)exp _ 

TT V ^0 



a3 



Tc(BEC) = 2 



(/^)|co| 

(44) 

r(3/2)C(3/2)J hJ' 



respectively [48l]. Here, 7 is Euler's constant and Dq{ij) is the density of states. In our low- 
density neutron matter, |co|/(a^t) is 5-7, and corresponds to the middle region in Fig. [151 
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FIG. 14: ^5*0 phase diagram of low-density neutron matter. The solid and open symbols with 
statistical uncertainties show the LO and NLO results, respectively. The dotted curves for Tc 
and T* are drawn by extrapolation. Neutron matter is in the superfluid phase below the critical 
temperature of the second-order phase transition. Above Tc, neutron matter is in the pseudogap 
phase 



49l |. in which pairing remains locally without forming long-range order, and undergoes a 
smooth transition from the pseudogap phase to the normal phase around T*, as pairing gets much 
less. 



The figure clearly shows that the thermal property of low- density neutron matter is not in a 
state of BCS, hut of BCS-BEC crossover. Though not discussed here, the Cq dependence of 



T* also verifies this point 471. |48|. 



The preceding point is perhaps better clarified by the cq dependence of the chemical 
potential /i. /x is positive in the weak-coupling BCS region and becomes negative in the 
strong-coupling BEC region by exhibiting a bosonic nature. Figure [16] illustrates the cq 
dependence of /i in the LO calculation, decreases as cq increases, and it takes a relatively 
small, positive value in the region of our low-density neutron matter. The small positive 
value is in accord with the neutron matter being close but not (yet) in the BEC region 
and indeed confirms the simple characterization of the crossover, a negative and small (in 
magnitude) value of l/ikpaQ) [28|, as noted in Sec. I. 
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FIG. 15: EFT parameter (cq) dependence of the critical temperature Tc- For easier comparison, 
Tc and cq are expressed as dimensionless by use of the spatial lattice spacing a and the hopping 
parameter t. The open circles are shown for Ns = 6^ at the quarter-filling (n = 0.5). The dashed 
curves are Tc/t at the BCS and BEC limits of Eq. 
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FIG. 16: Chemical potential as a function of the interaction strength cq in a dimensionless unit, 
with the spatial lattice spacing o and the hopping amplitude t. The calculation is of the LO for 
= 6^ and n = 0.5. 
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B. Pairing gap A 




0.2 0.4 0.6 0.8 1 

kp (fm-1) 

FIG. 17: Comparison of our Monte Carlo A to other calculations as a function of the neutron 
matter density (represented by the Fermi momentum kp). The solid diamonds show our results, 
with statistical uncertainties. The other calculations consist of three types: quantum Monte Carlo 
(symbols with statistical uncertainties), BCS (solid curve), and BCS with higher-order effects (R's, 
C's, and RG; shown by dotted and dashed curves). See text for the description of each calculation 
shown. 



Figure [T7] illustrates the density dependence of various A's reported in the literature. 
A's in the figure consist of those obtained by three types of calculations: (1) BCS (shown 
by a solid curve), (2) BCS or similar approximations, with higher order effects (dotted and 
dashed curves), and (3) quantum Monte Carlo (shown with error bar symbols). 

(1) Below kp ~ 0.7 fm~^, there are few recognizable differences 10, lUl among Abcs's 
calculated by various conventional A^A^ potentials: Argonne vis STj, Nijmegen 58|, and CD 



Bonn 



59| . Accordingly, Abcs's are represented by a single (solid) curve in Fig. [T71 



(2) Figure [T7] includes A's by the calculations beyond BCS. Calculations in the random 
phase approximation (RPA) with polarization effects are by Wambach et al. jssi (denoted as 
Rl), by Schulze et al. G^] (R2), and by Cao et al. 61] (R3). Calculations using correlated- 
basis functions are by Chen et al. (CI) and by Fabrocini et a l. (l6 | (C2). A calculation 
based on a renormalization group approach is by Schwenk et al. [631] (RG). The curves for 
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these A's are taken from similar figures in the recent hterature: Figs. 1 and 2 of Ref. 17 1 



and Fig. 4 of Ref. 



141 ]. In addition, though not shown, an extrapolation from finite nuclei 



results obtained by Hartree-Fock-Bogoliubov calculations also gives A close to the Abcs for 



kp < 0.5 fm 



-1 



641 ]. We see that these A's differ appreciably among each other, though 



recent works tend to give the values closer to the BCS A. 



(3) Two 



GFMC 



13 



:ypes of quantum Monte Carlo calculations have been reported based on the 



14l ] and AFDMC [l3] methods. The two methods are applied for a fixed 
number of neutrons using the conventional A^A^ potentials (or some model potentials), while 
our work is based on a grand canonical ensemble formulation. Figure [17] shows the most 
recent results of the GFMC [l^ (open squares), the AFDMC 17] (open circles), and ours 
(taken from Table HVl and shown by solid diamonds). 

In the figure, we see that all quantum Monte Carlo calculations are, overall, close to the 
Abcs- The AFDMC A is quite close to the Abcs in the density region examined in this 
work, while the GFMC A is smaller than the Abcs and is similar to (even slightly lower 
than) our A. Note that above kp ^ 0.6 fm~^, the AFDMC A becomes quickly smaller than 
the Abcs as the density increases. 

It is difficult to assess the three quantum Monte Carlo calculations by comparing them 
because the intermediate steps of the calculations are all different. Here, however, we point 
out a possible issue closely tied to their basic formulations and setups: stemming from the 
neutron numbers being fixed, the GFMC and AFDMC A's are calculated using the odd-even 
staggering (or the second-order finite difference) of the energy per neutron, 

A(odd N) = E{N) - ^ [E{N - 1) + E{N + 1)] , (45) 

where A^ is the number of neutrons. As described in Sec. Ill, our A's are calculated directly 
from the spin pair-pair correlation functions. By physical arguments, the two ways of cal- 
culating A are expected to be the same for a large A^, but we are not aware of a rigorous 
proof for this expectation. Since it has been a common practice to apply Eq. fj45l) for the 
extraction of A from finite nuclei [l|, [s^] , closer examination of this issue would be desirable. 



as exemplified in Ref. 64]. 



As noted above, it is desirable to apply Eq. 



A^ = 92 are used in the GFMC calculation 



5]) for a large A^. The large values up to 



J4, while up to A^ = 68 in the AFDMC 



17|. 



Both A^'s are perhaps large enough to provide reliable information for A^ — > oo. While it 
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might be caused by the different ways the nuclear potentials are applied in the two methods, 
the noticeable difference between the GFMC and AFDMC A's is puzzling to us. 



C. Further improvement of the present work 

We note here the aspects of this work that we would like to improve. 

(1) The largest lattice size we have used is A^^ = 10^, but larger lattices would be desirable 
for reliably reaching the thermodynamic limit. For this, we would like to study more closely 
the use of the hybrid Monte Carlo (HMC) method. As the commonly used method in 



lattice QCD calculations 



the HMC is expected to reduce the computation time from 



(NsNt)^ or {NJVtY (for the DQMC) to ~ {NsNtf/^. Our trial apphcation of the HMC 



(following Ref. 65|]) in our problem has shown a strong dependence on the HMC parameters, 
such as the size and number of molecular dynamics steps and has brought about a difficult 
compromise between the computation time and the systematic error. We suspect that the 
difficulty stems from badly conditioned fermion matrices and also from our (effectively) 
strong interaction. We would like to resolve this issue and find a practical procedure for 
optimizing the HMC calculation for this problem. 

(2) Because of lack of time, we have examined the continuum limit by applying the case 
of Ng = 6^ to all A^s's that we computed. The possible A^, dependence is a potentially 
important source of the systematic error, and we would like to clarify this issue. 

(3) The matching of the LO and NLO calculations indicates that our A deviates from the 
Abcs more appreciably in the matching density region, kp ~ 0.15-0.3 fm~^. It is difficult to 
establish the deviation by using the present statistics. We would like to examine this density 
region more closely to determine whether such a fine structure of the density dependence of 
A exists. 



VIII. SUMMARY 

In conclusion, we have investigated thermal properties of low-density neutron matter 
by the determinantal quantum Monte Carlo lattice calculations with the single- and two- 
parameter pionless EFT A^A^ potential. The ^5*0 pairing gap at T 0, the critical tem- 
perature of normal-to-superfiuid phase transition, and the pairing temperature scale have 
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been determined directly from the correlation functions and have provided the temperature- 
density phase diagram for the density of (10~^-10~^)po- The thermodynamic hmit was taken, 
and the continuum limit was examined in the determination. The pairing gap was found 

to be approximately 30% less than the BCS value. The physics of neutron matter in this 
density region has clearly been identified as a BCS-BEC crossover. 
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APPENDIX A: DETERMINATION OF THE EFT POTENTIAL PARAMETERS 

Co (A) AND C2(A) 



The EFT potential parameters, Cq{A) and C2(A), are determined from the observables for 
an appropriately chosen value of A. As the observables, we choose the scattering length 
and the effective range tq in the effective range expansion of Eq. ([6]) with A = n/a in our 
lattice calculation (where a is the lattice spacing). 

A is needed in the determination of Co(A) and C2(A) so as to regularize loop contributions, 
which otherwise diverge. With the regularization, the Schrodinger equation is solved, and 



and To are expressed in terms of Co(A) and C2(A) algebraically [37|, l66|. The direct 
use of the algebraic expressions, however, amounts to a mere phenomenological fit. As an 
application of EFT, we must ensure that EFT counting rules are properly applied: because 
our EFT Lagrangian is truncated at p^/Q"^, we must be consistent with the truncation in 
the determination of Cq and C2. That is, C2(A) must be treated perturbatively by neglecting 



the 0([c2(A)]^)-order contributions. We then obtain [371] 



M 1 
4:71 ao 



1 M , 
+ —Li 



Co (A) 27r 



M C2(A) 

TC^ Co(A) 



M C2(A) M 1 , , , , 



167r " cl{A) Att^A 
where Li = 9iA and L3 = 9^A^. The numerical values of 61,63, and -R(O) for large lattices 
are given in Ref. js^]- The inversion of Eq. (]A1|) is, again by treating C2(A) perturbatively. 



co(A) = cr(A)|l + ^(|)'L3,[cr(A)]n-cr(A) + Aco(A), 

C2(A) = ^v[cPiA)r, (A2) 
where r] = 1 + 4R{0)/ (vrroA), and the leading-order co(A), c^^\A), is given by 

4"W-'^(---Ly. (A3) 



M \ao 71 

Equations (lAip and (lA2p consistently include up to their combined use is 

equivalent to solving the Schrodinger equation with the truncated potential of Eq. flA2p 
by treating C2(A) perturbatively. That is, in this treatment, we obtain exactly the same 
Oo and ro as those determined phenomenologically or obtained by solving the Schrodinger 
equation with no counting rule applied. Because of this, the phase shifts determined by 
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Oo and tq are also exactly the same as those determined by the LO and NLO potentials 
by consistently applying the EFT counting rule. The same EFT treatment should also be 
applied to calculations of many-nucleon systems, as we have done in this work. Note that 
upon the application of the EFT counting rule, consistency is the vital point, as is evident 



from the observation 
not properly applied 



;hat tq turns out to be negative for a certain range of A if this step is 



66]. 
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68|. 
36|. 



For ao and tq, we have used the old values of —16.45 fm and 2.83 fm, respectively 
The most recent values are Oq = — 18.9±0.4 fm and tq = 2.75±0.11 fm as quoted in Ref. 
The discrepancy between the two ao values is 13 ± 2 % and not negligible, but its effects are 
expected to be much smaller. 

As Eq. f lA3p implies, c^^ is dominated by the A contribution because c^q^ is close to 



3G 



691], dictated by the large 



the nontrivial fixed point in the renormalization group flow 
magnitude of Oq. Consequently, Cq and C2 are quite insensitive to the exact value of Oq. For 
example, at kp = 60 MeV, using the standard parameter set of Table I, the NLO co/(a^t) 
and C2/{a^t) differ by 1.8% and 1.4%, respectively, between the old and the most recent 
values of ao and tq. The corresponding LO co/(a'^t) differs by 1.4% between them. 

Generally some A contributions must cancel in calculating observables, so that their 
values are independent of the regularization procedure. But the closeness to the fixed point 
suggests the cancellation to be effectively small in this case. Although repeating our entire 
calculations is quite time consuming and unrealistic at present, we have performed a limited, 
test LO calculation at kp = 60 MeV for A^^ = 6^ and n = 1/4. We find A differs by about 2%, 
in the same order of the statistical uncertainties of the Monte Carlo calculation: A = 0.63(1) 
and = 0.64(3) MeV for = —16.45 and = —18.9 fm, respectively. This finding also confirms 



the following observation: in the accompanying paper [29|], we report the determination of 
various quantities at the unitary limit (|ao| oo with tq = 0) by making the extrapolation 
7] = l/laokp) — > 0. By taking the rj variation to be an ao variation, we find that the above 
discrepancy in A is 2.2% for kp = 60 MeV and decreases as kp gets larger and increases as 
kp gets smaller. 
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APPENDIX B: PHYSICAL SIZES OF A NEUTRON PAIR AND COMPUTA- 
TIONAL LATTICE 



A measure of the size of an interacting neutron pair (a Cooper pair) in the superfluid 



state J ^cp; is 



271 



= (Bl) 

^cp niust be smaller than the dimension of the cubic lattice, as a necessary condition for the 
simulation of the collective state (but clearly not a sufficient one). Table |VT] shows that ^cp 

1/3 

is indeed much smaller than the dimension of the lattice aNs , except for the marginal case 
of Ng = 4^. Note that ^cp depends on a and L through the n dependence of A. The a and 
Ng dependence of ,^cp through A is weak, as seen in Sec. IV A. In the table, we list ,^cp for 
Ns = 4^ and n = 1/4, for simplicity. 

TABLE VI: Physical sizes of a neutron pair and computational lattices. 



kp (MeV) Cap (fm) a (fm) aN^'^ {Ns = 4^) aN}'-^ (iV, = 6^) aN}'-^ {Ns = 8^) aiYs^ {Ns = 10^) 



15 


1.3 X 10^ 


25.64 


102.6 


153.8 


205.1 


256.4 


30 


47 


12.82 


51.3 


76.9 


102.6 


128.2 


45 


28 


8.55 


34.2 


51.3 


68.4 


85.5 


60 


21 


6.41 


25.7 


38.5 


51.3 


64.1 


90 


11 


4.27 


17.1 


25.6 


34.2 


42.7 


120 


8.6 


3.21 


12.8 


19.2 


25.6 


32.1 



APPENDIX C: TECHNICAL DETAILS OF MONTE CARLO COMPUTATION 

In this appendix, we discuss some technical details of the setup for the implementation 
of our lattice calculations. 
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1. Parameter values 



The parameter set for lattice sizes is the following: the number of spatial lattice sites 
used are 

Ns = 4^ 6^ 8^ and 10^ (CI) 

so as to extrapolate the data into the thermodynamic limit (A^^ oo); the number of 
temporal lattice sites is 

4<Nt< 128, (C2) 



where the discretization size of the temporal lattice is the same in Ref. 47| as 



AN, . (C3) 

The typical example of one production run is as follows. Because the method of grand 
canonical ensemble is used, fi is fixed in each run. The thermal observable for the desired 
density p is interpolated from a few sets of the observables calculated at different fi. About 
1000-10000 samples are accumulated to obtain statistics with a precision of several percent. 

2. Determinantal quantum Monte Carlo 

a. Temporal lattice spacing 

To choose A/5, we need to know how the expectation values of thermal observables are 
affected by the choice. Figure [T8] illustrates the dependence of Af3 on the thermal observable 
Ca in our DQMC calculation. The data have been taken with fi/t = and A^^ = 4^ at 
kp = 30 MeV. The figure is a typical example, and we have observed similar results with 
other thermal observables and parameter values. 

From Fig. [181 we see that the expectation values of thermal observables are affected only a 
little for A/3 ^ 0.2t, confirming that the choice employed in the previous DQMC calculation 



similar to ours 



201 ] is indeed reasonable, and so we adopted this choice. 



b. Prethermalization steps 

At the start of sampling, we generate the initial configuration of the auxiliary fields x- Iii 
our DQMC calculation, we use the hot start, in which a random (disordered) configuration 
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FIG. 18: Pair correlation function Ca as a function of temporal lattice spacing A/3, with Ng = 4^ 
at kp = 30 MeV in the unit of hopping amplitude t in the DQMC calculation. 

is used, instead of the cold start using a uniform (ordered) configuration. Following the 
start, we must take a sufficient number of prethermalization steps to obtain the equilibrium 
configurations, statistically independent from the initial configuration in the Markov chain. 

Figure [19] illustrates the dependence of the sample number on the thermal observable Ca 
in our DQMC calculation. The data have been taken with fi/t = —1.83 at Ng = and 
Nt = 12. The figure shows that the equilibrium starts to be reached after 100-150 samples. 
Similar results are observed with other observables and for other parameter values. 



c. Thermalization steps and autocorrelations 

To ensure statistically independent configurations, we must take thermalization (decorre- 
lation) steps between sample takings. We determine the number of the thermalization steps 
by monitoring the autocorrelation. The autocorrelation for k conservative samples of the 
observable O, Co{k), is of the standard form 



Co{k) 



(C4) 



where (■ ■ ■ ) denotes the average over the random walk labeled with i, for example, 

^ N-k 

(OiO.+k) = 0{X,)0{X,^,). (C5) 



i=l 
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FIG. 19: The pair correlation function Ca as a function of the sample number at A'^s = 4^ and 
Nt = 12 in the unit of hopping amplitude t in our DQMC calculation. 

The condition of no correlation is Cq ~ 0, but in practical terms Cq ^ 0.1 is recommended 
[3], and thus we ensure Co to be less than 10% . 

A typical case of the autocorrelations for some observables is shown in Fig. [20] with the 
parameter set {Ng = 4^, Nt = 12, and kp = 30 MeV). The autocorrelations are seen to be 
less than 0.1 for more than ten thermalization steps between samples. 

3. Systematic error of the DQMC 

Here, we discuss the systematic uncertainties of the DQMC besides the statistical ones 
due to data sampling. After ensuring the independence between samples by keeping the 
autocorrelations of thermal observables small enough as described in Appendix A 2, the 
systematic error of the DQMC on observables solely comes from the size of the discretization 
of the time slice Af3, which is related to the inverse of temperature (3 = NtAf3. 

For confirming the consistency of our DQMC calculation with others, we compare Tc/t 



with that in Refs. 



47 



48l | over the various interaction strengths Co/(ci^t) at fixed temporal 



lattice spacing A/5 = 0.125/t, which has been commonly used in the condensed-matter 
physics. For estimating the systematic errors caused by finite A/3, the A/3 dependence of 



A, Tc, and T* have also been further examined. 
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FIG. 20: Autocorrelation as a function of thermalization steps between samples taken with the 
number of spatial lattice sites A'^, = 4^ and of temporal lattice sites Nt = 12 at the Fermi momentum 
kp = 30 MeV in our DQMC calculation. 

By these preliminary DQMC calculations, we can ensure the consistencies of DQMC 
calculations with those in other literature. The systematic uncertainties caused by our 
calculations with finite A(3 amount to around 10%. 



a. Comparison of Tc{co/{a^t)) with other work 

First we ensure that our DQMC calculation at finite Af3 is consistent with other literature. 
Figure [T5] is the critical temperature Tc as a function of interaction strength \co\/{aH) at the 
quarter-filling (n = 1/2) in A^^ = 6'^. is obtained through the infiexion point of the curve 
of pair correlation function Ca- The parameters used in the calculations are Af3 = 0.125/t, 



-^pretherm — 200, A< 



therm 



50, A'sampic = 1000-2000. Our T^{\co\/ (aH)) over the interaction 



strength ranging between BCS and EEC limits is in good agreement with Refs. 47|, |48|] of 
the same setup within around 5% of errors, which is within the DQMC results in other 
literature, ranging around 10% at half-filling {n = 1) as shown in the left panel of Fig. 5.13 



in Ref . 



48|. 
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b. Dependence of thermal observables on A/3 



Now that our DQMC calculations with finite A/3 are confirmed within around 5% of the 
differences, we have to consider the systematic error from the discretization of temporal 
direction A(3. Figure shows the dependence of various thermal observables on A(3 by 
fixing T/t = l/{NtAf3t) = 0.4. The expectation values of thermal observables are obtained 
by 1000-2000 samples with A^pretherm = 200 and A^therm = 100 at the one-eighth filling 
(n = 1/4). In Fig. [21], we take the ratio of thermal observables at Af3 = 0.125/t to those 
at the continuum limit of the temporal direction A/5 ^ to make the deviations easily 
visible. As summarized in Table IVTII the differences of the observables with At = 0.125/t 
and A/3 are around 5% (for xp), 10% (for Ca and E/A), and 20% (for /x). Note that 
we use only Ca and xp ioi obtaining Tc and T* in this work. 
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FIG. 21: A/3 dependence of the ratio of energy per particle E/A, pair correlation function Ca, 
Pauli spin susceptibility xp, and chemical potential /j, to those at A/3 — > at T/t = 0.4 with the 
interaction streng th co/{a^t) = -6.0 at the one-eighth filling (n = 1/4) in the dimensionless unit. 



Next we examine the influence of flnite A/3 on and T*. Figures [221 and [23] summarize 
the effect of the flnite A(3 on and T*. As seen in those flgures, T^/t = 0.45(1) MeV and 
T*/t = 0.87(2) MeV for A/3 = 0.125/t, and Tjt = 0.47(2) MeV and T*/t = 0.79(2) MeV for 
A/3 = 0.0625/t. The quantities in the parentheses indicate the statistical uncertainties. The 
deviations in and T* without the statistical errors are around 5% and 10%, respectively. 
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TABLE VII: Ratio of thermal observables. 



O 0{A(3 = 0.125/t) 0{Ap 0) 0{A/3 = 0.125/i)/O(A/3 0) 



E/{At) 


1.449(8) 


1.625(8) 


0.892(9) 


Ca 


1.45(1) 


1.32(1) 


1.10(2) 


XP 


0.207(9) 


0.1965(9) 


1.05(5) 


/x/t 


1.49(1) 


1.23(1) 


1.21(1) 



We have to count on these discrepancies of around 10% as the systematic error of our final 
results besides the statistical error. 



6 
5 
4 



< or 

O 3 r 



2 r 
1 - 



—I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 



S 



i 



— I — I — I — I — I — I — I — I — 

ANt = 0.125/t 



AN, = 0.0625/t 



J , , , , L 



0.5 



1 

T/t 



J , , , L. 



1.5 



FIG. 22: Pair correlation function Ca as a function of temperature T in the unit of hopping 
amplitude t at different A/? at kp = 30 MeV, iVg = 4^, and n = 1/4. The open and solid circles 
with statistical errors are the results at LO and NLO, respectively. 



As described in Sec. Ill, we use Pg for an estimation of A. The constant tails of Pg at 
the large separation of pairs are Ps(A/3 = 0.125/t) = 0.02784(46) and Ps(A/3 = 0.0625/t) = 
0.0295(24) at kp = 30 MeV. The resultant pairing gaps extracted from Pg through A = 
cov^ with Co = 0.8012 MeV are A(A/3 = 0.125/t) = 0.1337(11) MeV and A{A(3 = 
0.0625/t) = 0.1377(56) MeV. The deviation between them without the statistical errors 
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FIG. 23: Pauli spin susceptibility xp ^ ^ function of temperature T in the unit of hopping 
amplitude t at different A(5 at kp = 30 MeV, Ns = 4^, and n = 1/4. The open and solid circles 
with statistical errors are the results at LO and NLO, respectively. 

quoted by the parentheses is 0.004 MeV, which results in around 3% of the systematic error. 
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